Asymmetric Abelian Avalanches and Sandpiles

We consider two classes of threshold failure models, Abelian avalanches and sandpiles, with the redistribution matrices satisfying natural conditions guaranteeing absence of infinite avalanches. We investigate combinatorial structure of the set of recurrent configurations for these models and the corresponding statistical properties of the distribution of avalanches. We introduce reduction operator for redistribution matrices and show that the dynamics of a model with a non-reduced matrix is completely determined by the dynamics of the corresponding model with a reduced matrix. Finally, we show that the stationary distributions of avalanches in the two classes of models: discrete, stochastic Abelian sandpiles and continuous, deterministic Abelian avalanches, are identical. Introduction. Different cellular automaton models of failure (sandpiles, avalanches, forest fires, etc.), starting with Bak, Tang and Wiesenfeld [BTW1, BTW2] sandpile model, were introduced in connection with the concept of self-organized criticality. Traditionally, all of these models are considered on uniform cubic lattices of different dimensions. Recently Dhar [D1] suggested a generalization of the sandpile model with a general (modulo some natural sign restrictions) integer matrix ∆ of redistribution of accumulated particles during an avalanche. An important property of this Abelian sandpile (ASP) model is the presence of an Abelian group governing its dynamics. Dhar introduced the set of recurrent configurations for an Abelian sandpile model, the principal geometric object governing its dynamics in the stationary state. The burning algorithm introduced in [MD] allows to recognize, for a symmetric sandpile model, when † Currently at Department of Mathematics, University of Toronto, Toronto M5S 1A1, Canada Asymmetric Abelian Avalanches and Sandpiles Page 2 a stable configuration is recurrent. A more sophisticated script algorithm suggested in [S] plays the same role for asymmetric models. Both algorithms provide, in fact, certain information on the combinatorial structure of the set of recurrent configurations. Abelian sandpiles were also studied in [WTM, D2, DM, GM, Cr, Go, CC, G1, G2, S, DMan]. In a non-dissipative case ( ∑ j ∆ij = 0, for all i) an avalanche in the ASP model coincides with a chip-firing game on a directed graph [BLS, BL] where −∆ is the Laplace matrix of the underlying digraph. Another class of lattice models of failure, slider block models introduced by Burridge and Knopoff [BK] and studied in [CaL, Ca, N, MT, LKM], as well as models [FF, D-G, OFC, CO, Z, PTZ, GNK] which are equivalent to quasistatic block models, have continuous time and continuous quantity at the lattice sites which accumulates in time and is redistributed during avalanches. This quantity is called the slope, height, stress or energy by different authors. In slider block models it corresponds to force [OFC]. We use the term height as in [D1]. Gabrielov [G1], introduced Abelian avalanche (AA) models, deterministic lattice models with continuous time and height values at the sites of the lattice, and with an arbitrary redistribution matrix. For a symmetric matrix, these models are equivalent to arbitrarily interconnected slider block systems. In the case of a uniform lattice, these models were studied in [FF, D-G] and in [GNK] (as series case a). The stationary behavior of the AA model is periodic or quasiperiodic, depending on the loading rate vector. At the same time, the distribution of avalanches for a discrete, stochastic ASP model is identical to the distribution of avalanches for an arbitrary quasiperiodic trajectory (or to its average over all periodic trajectories) of a continuous, deterministic AA model with the same redistribution matrix and loading rate [G1]. In this paper, we introduce general conditions on redistribution matrices that are equivalent to the absence of infinite avalanches in the model. The models satisfying these conditions include, in particular, the models with non-negative dissipation or codissipation Asymmetric Abelian Avalanches and Sandpiles Page 3 considered before. We allow spatially inhomogeneous loading rate and show that the set of recurrent configurations does not depend on the loading rate vector, as long as a natural condition guaranteeing absence of non-loaded components in the model is satisfied. We continue the study of the combinatorial structure of the set of recurrent configurations started in [D1, DM, MD, S] and introduce a new description of this set, essentially improving the script algorithm suggested in [S] for asymmetric redistribution matrices, Additional possibility for the study of the models with asymmetric matrices arises from the matrix reduction operations. These operations act on the redistribution matrices in the same way as the topple operations act on unstable configurations, and satisfy the same property of the independence of the resulting reduced matrix on the possible change of the order of reductions. Each reduction operation simplifies the redistribution matrix, and replaces the original model by a simpler reduced model such that the combinatorics of the set of recurrent configurations for the reduced model completely determines the combinatorics for the original non-reduced model. In the first section, we define configurations, redistribution matrices and avalanche operators. We show, following a construction implicitly present in [S], that a legal sequence of topples satisfies certain minimality condition among all (possibly illegal) sequences of topples with the same final stable configuration (lemma 1.1). This minimality condition provides, in particular, a new proof of the principal Abelian property (theorem 1.2). We introduce the class of avalanche-finite redistribution matrices satisfying eight equivalent conditions and check these conditions for the matrices with non-negative dissipations and codissipations. In the second section, we define, following [G1], the AA model as a sequence of loading periods and avalanches and describe the dynamics of the model on its attracting set of recurrent configurations. The arguments here are similar to the arguments of Dhar [D1] for the ASP models. In the third section, we study the combinatorial structure of the set of recurrent Asymmetric Abelian Avalanches and Sandpiles Page 4 configurations of the AA model. The principal result here, the theorem 3.8, describes this set as the complement in the set of all stable configurations to the union of negative octants with the vertices in a finite set N . An explicit constructive description of this set N is given in the theorem 3.11. In the fourth section, we show the possibilities to extract the information on the dynamics of an AA model from the dynamics of another model with a simplified redistribution matrix. We introduce the total reduction operator for redistribution matrices, similar to the avalanche operator for configurations. We show that the stationary dynamics of an AA model with a redistribution matrix ∆ is completely determined by the dynamics of the corresponding model with a reduced redistribution matrix, the total reduction of ∆. In the fifth section, we introduce marginally stable configurations and derive formulas for the mean number of avalanches. The arguments here are again similar to those of Dhar [D1], modified for the more general situation considered here. In the sixth section, we establish the identity between the distributions of avalanches for AA and ASP models with the same redistribution matrices. Some of the results of this paper were announced in [G2]. 1. Redistribution matrices and avalanches. Let V be a finite set of N elements (sites), and let ∆ be a N × N real matrix with indices in V . We call ∆ a redistribution matrix when ∆ii > 0, for all i; ∆ij ≤ 0, for all i 6= j. (1) A real vector h = {hi, i ∈ V } is called a configuration. The value hi is called the height at the site i. For every site i, a threshold Hi is defined, and a site i with hi < Hi is called stable. A configuration is stable when all the sites are stable. For i ∈ V , a topple operator Ti is defined as Ti(h) = h− δi (2) Asymmetric Abelian Avalanches and Sandpiles Page 5 where δi = (∆i1, . . . ,∆iN) is the i-th row vector of ∆. Obviously, every two topple operators commute. The topple Ti(h) is legal if hi ≥ Hi, i.e. if the site i is unstable. No topples are legal for a stable configuration. A sequence of consecutive legal topples is called an avalanche if it is either infinite or terminates at a stable configuration. In the latter case, the integer vector n = {ni, i ∈ V } where ni is the number of topples at a site i during the avalanche is called its script, and the total number ∑ i ni of topples in the avalanche is called its size. The following lemma shows that the avalanches are “extremal” among all the sequences of (possibly, illegal) consecutive topples with the same endpoints. It allows, in particular, to give an alternative proof of the principal property of avalanches — the script and the final stable configuration depend only on the starting configuration, not on the possible choice in the sequence of topples (theorem 1.2 below). Lemma 1.1. Let h be an arbitrary configuration and let m be an integer vector with non-negative components mi such that g = h− ∑ miδi is a stable configuration. For any finite sequence of consecutive legal topples started at h, with ni topples at a site i, we have mi ≥ ni. Proof. The arguments appear implicitly in [S]. We use induction on the size n = ∑ ni of the sequence of legal topples. For n = 0, the statement is trivial. Let it be true, i.e. mi ≥ ni, for a sequence with ni topples at a site i. If a site j is unstable for a configuration f = h − ∑ niδi then gj < fj . Due to (1), this implies mj > nj, hence the statement remains true when we add a topple at the site j to the sequence. Theorem 1.2. (Sf. [D1], [BLS], [BL].) Every two avalanches starting at the same configuration h are either both infinite or both finite. In the latter case, the scripts of both avalanches coincide. In particular, both avalanches terminate at the same stable configuration and have the same size. Proof. The statement follows easily from the lemma 1.1. Asymmetric Abelian Avalanches and Sandpiles Page 6 Remark 1.3. If we consider a configuration as an initial state of a game, and every legal topple as a legal move, an avalanche becomes a (solitary) game. The theorem 1.2 means that this game is strongly convergent in the definition of [E]. Lemma 1.4. For every site i that toppled at least once during an avalanche, hi ≥ Hi−∆ii till the end of the avalanche. The statement follows from (1) and (2). Let R+ = {hi ≥ 0, for all i} and R− = −R+ denote positive and negative closed octants in R , and let Ṙ+ = {hi > 0, for all i} be an open positive octant. Let ∆′ be the transpose of the matrix ∆. Theorem 1.5. For a redistribution matrix ∆, the following properties are equivalent. i. Every avalanche for ∆ is finite. ii. ∆(R+ \ {0}) ∩R− = ∅. iii. ∆(R+) ⊇ R+, i.e. ∆−1 exists and all its elements are non-negative. iv. ∆(R+) ∩ Ṙ+ 6= ∅. i′. Every avalanche for ∆′ is finite. ii′. ∆(R+ \ {0}) ∩R− = ∅. iii′. ∆(R+) ⊇ R+, i.e. ∆′−1 exists and all its elements are non-negative. iv′. ∆(R+) ∩ Ṙ+ 6= ∅. Proof. (ii′)⇒ (i). Let us show that for ∆ satisfying (ii′), every avalanche is finite. If there exists an infinite avalanche started at a configuration h, let r(k) = {ki/k, i ∈ V } where ki is the number of topples at a site i after a total number of topples k. According to (2), the configuration after k topples is h(k) = h − k∆′r(k). Let r ∈ R \ {0} be an accumulation point for r(k) (it exists because all these vectors have unit length) and p = −∆′r an accumulation point for (h(k)− h)/k. According to lemma 1.4, components of h(k)− h are bounded from below. Hence all components of p are non-negative, and ∆ does not satisfy (ii′). Asymmetric Abelian Avalanches and Sandpiles Page 7 (i)⇒ (iii′). Let h be a configuration in R+, and let a finite avalanche starting at kh terminates at a stable configuration h(k). Let r(k) = {ki/k, i ∈ V } where ki is the number of topples at a site i during this avalanche. We have ∆′r(k) = h− h(k)/k. Let r be an accumulation point for r(k), as k → ∞. Then r ∈ R+ and ∆′r = h, because h(k) remains bounded as k→∞. (iv)⇒ (ii′). Suppose that ∆ does not satisfy (ii′). This means that there exists a linear form l 6= 0 with non-negative coefficients such that l(δi) ≤ 0, for all i. Hence l( ∑ ciδi) ≤ 0 for any combination of the vectors δi with non-negative coefficients ci. At the same time, l(δ) > 0, for every δ ∈ Ṙ+. This means that ∆(R+)∩ Ṙ+ = ∅ and ∆ does not satisfy (iv). (iii′)⇒ (iv′). The implication is obvious. Combining the four implications, we have (iv)⇒ (ii′)⇒ (i)⇒ (iii′)⇒ (iv′). The same arguments applied to ∆′ instead of ∆ imply (iv′)⇒ (ii)⇒ (i′)⇒ (iii)⇒ (iv). This completes the proof. Definition 1.6. A redistribution matrix satisfying the conditions of the theorem 1.5 is called avalanche-finite. Remark 1.7. Let ∆ be an avalanche-finite matrix, and let t ∈ Ṙ+, ∆t ∈ Ṙ+. Such a vector t always exists due to the property (iii) or (iv). Let |h|t = (h, t) be the t-weighted length of a configuration h. Then |Ti(h)|t < |h|t, for all i ∈ V , i.e. every topple operator dissipates the t-weighted length. This can be also used to prove the implications (iii)⇒ (i) and (iv)⇒ (i). Definition 1.8. The value si = ∑ j ∆ij is called the dissipation at the site i, and the value sj = ∑ i ∆ij is called the codissipation at the site j. A site i is called dissipative (non-dissipative) if si > 0 (si = 0). A site j is called codissipative (non-codissipative) if sj > 0 (s ′ j = 0). An underlying digraph Γ = Γ(∆) of a redistribution matrix ∆ is defined by the vertex set V (Γ) = V and an edge from a site i to a site j drawn iff ∆ij < 0. Asymmetric Abelian Avalanches and Sandpiles Page 8 Let s′ be a diagonal matrix with sii = s ′ i, and let ∆0 = ∆− s′ (3) be the non-codissipative part of ∆. The matrix ∆0 coincides with the Kirchhoff matrix of Γ, with conductance of an edge −→ ij defined as −∆ij [T, p.138]. A subset W of V is called a sink in Γ if there are no edges from sites in W to sites outside W , and a source if there are no edges from sites outside W to sites in W . A matrix ∆ is called weakly dissipative if all the dissipation values si are non-negative and the digraph Γ(∆) has no non-dissipative sinks, i.e. from every site there exists a directed path in Γ(∆) to a dissipative site. Proposition 1.9. A matrix with non-negative dissipation values is avalanche-finite if and only if it is weakly dissipative. Proof. If the graph Γ(∆) has a non-dissipative sink W ⊆ V then ∑ i∈W hi does not decrease during an avalanche, hence the avalanche started at a configuration with large enough values of hi, i ∈W , cannot be finite. Suppose now that ∆ is weakly dissipative. It follows from the definition 1.8 that σ = ∑ i hi does not increase at any topple and decreases when a dissipative site topples. Suppose that there exists an infinite avalanche, and let W ⊂ V be the subset of sites that topple infinite number of times in this avalanche. Then all the sites in W are nondissipative, otherwise σ would decrease indefinitely, in contradiction to the lemma 1.4. At the same time, W is a sink of Γ, otherwise hj would increase indefinitely at any site j 6∈W such that ∆ij < 0, for some i ∈W . This contradicts the definition 1.8. Definition 1.10. A matrix ∆ is called weakly codissipative if all the codissipation values sj are non-negative and the digraph Γ(∆) has no non-codissipative sources, i.e. to every site there exists a directed path in Γ(∆) from a codissipative site. Proposition 1.11. A matrix with non-negative codissipation values is avalanche-finite if and only if it is weakly codissipative. Asymmetric Abelian Avalanches and Sandpiles Page 9 Proof. The statement follows from the theorem 1.5 and proposition 1.9, because the transpose of a weakly codissipative matrix is weakly dissipative. Proposition 1.12. For every avalanche-finite matrix ∆, det(∆) > 0. Proof. Let a redistribution matrix ∆ satisfy the condition (ii) of the theorem 1.5. For t ∈ [0, 1], all the matrices ∆t = t∆ + (1− t)E from the segment connecting with the unit matrix E satisfy (ii). Hence all these matrices are avalanche-finite. Due to the condition (iii), all the matrices in this segment are non-singular, hence their determinants have the same (positive) sign. In the following, we consider only avalanche-finite redistribution matrices. Definition 1.13. For a configuration h, the avalanche operator Ah is defined as the stable configuration that terminates an avalanche initiated at h. Due to the theorem 1.2, this stable configuration is unique. If h is stable, Ah = h. Example 1.14. The sandpile model introduced in [BTW], n× n square lattice with the nearest neighbor interaction and particles dropping from the boundary, is defined by a symmetric redistribution matrix ∆ of the size n × n. The rows and columns of ∆ are specified by a vector index i = (i1, i2) with 1 ≤ iν ≤ N , for ν = 1, 2, ∆i,i = 4, ∆i,j = −1, for i1 = j1, i2 = j2±1, and for i1 = j1±1, i2 = j2, ∆i,j = 0 otherwise. This matrix is weakly (co-) dissipative, hence avalanche-finite. Example 1.15. The 1-dimensional model with the failure depending on the local slope, introduced in [BTW] (for m = 1) and studied in [KNWZ, LLT, LT, S, CFKKP], is defined as follows. At every site i, 1 ≤ i ≤ N , we place ki particles, and set kN+1 = 0. The site i topples when ki − ki+1 ≥ m. The topple operator removes m particles from the site i and adds one particle to each site j = i+1, . . . , i+m as soon as j ≤ N . After the transformation hi = ki−ki+1, for 1 ≤ i ≤ N , this model can be defined by a redistribution matrix ∆ with ∆i,i = m + 1, for i < N, ∆N,N = m ∆i,i−1 = −m, for i > 1, ∆i,ν = −1, for i < N Asymmetric Abelian Avalanches and Sandpiles Page 10 and ν = max(i+m,N) (sf. [S]). This matrix is not symmetric when m > 1. It is weakly (co-) dissipative, hence avalanche-finite. Example 1.16. A chip-firing game introduced in [BLS, BL] is defined by a (directed) graph Γ with a certain number of chips placed at each of its vertices, and a sequence of legal moves (fires), when one particle is allowed to be moved from a vertex i to the end of each edge directed from i, in case the total number of chips at the vertex i is not less than the total number of the edges directed from i. The corresponding redistribution matrix is, after a sign change, the Laplace matrix of Γ. It is always degenerate (all the dissipations are equal 0) hence not avalanche-finite. Example 1.17. For N ≤ 3, a redistribution matrix ∆ is avalanche-finite iff det(∆) > 0. However, for N ≥ 4 there exist redistribution matrices with positive determinant which are not avalanche-finite. Consider, for example, an 4× 4-matrix ∆ =  1 −3 −1 0 −3 1 −1 0 −1 −1 1 −2 0 0 −2 1  . We have det(∆) = 16 > 0. At the same time, the matrix ∆ is not avalanche-finite, because ∆(1, 1, 0, 0) = (−2,−2,−2, 0), in contradiction to the condition (ii) of the theorem 1.5. 2. Abelian avalanche model. In this section, we define the Abelian avalanche model as a sequence of slow loading periods and fast redistribution events (avalanches). Many of the statements in this section are similar to the corresponding statements of Dhar [D1] for the ASP models. We present these statements with short proofs to make the paper self-contained. Also, the class of the redistribution matrices and loading vectors considered here is more general than in [D1]. Let v = {vi, i ∈ V } be a non-zero vector with non-negative components. For an (avalanche-finite) redistribution matrix ∆, an Abelian avalanche (AA) model [G1] with a loading rate vector v is defined as follows. Asymmetric Abelian Avalanches and Sandpiles Page 11 For every stable configuration, every height hi increases in time with the constant rate vi, until a height hi equals or exceeds a threshold value Hi at some site i. Then the site i topples according to (2) starting an avalanche which terminates at a stable configuration. After this the loading resumes, and the process continues indefinitely. Definition 2.1. An AA model is called properly loaded if the digraph Γ(∆) does not contain non-loaded sources, e.g. to every site there exists a directed path in Γ(∆) from a loaded site. Here a site i is called loaded if vi > 0 and non-loaded if vi = 0. If the model is not properly loaded, some parts of the system do not evolve in time. For a properly loaded model, it is easy to show that the rate of topples at every site is positive. The dynamics of the model does not change if we replace the values Hi by some other values, and add the difference to all configuration vectors. For convenience we take Hi = ∆ii. In this case, hi(t) ≥ 0 for any trajectory h(t) of the model when the i-th element has been toppled at least once. Hence, for a properly loaded model, only configurations h ∈ R+ are relevant for the long-term dynamics. Let S = {0 ≤ hi < ∆ii} be the set of all stable configurations in R+ . Remark 2.2. In the case of a symmetric matrix ∆ and vi = si, for all i, the AA model is equivalent to a system of blocks where i-th block is connected to j-th block by a coil spring of rigidity −∆ij , if ∆ij < 0, and to a slab moving with a unit rate by a leaf spring of rigidity si. For every block, a static friction force Hi is defined, and a block is allowed to move by one unit of space when the total force hi applied to this block from other blocks and the moving slab, equals or exceeds Hi. The weak dissipation, weak codissipation and proper loading conditions coincide in this case, and are satisfied when the system of blocks (including the moving slab) is connected. Proposition 2.3. (Sf. Dhar [D1].) Let r = {ri} be defined by ∆′r = v. Then ri is equal to the average, per unit time, rate of topples at a site i, independent of the initial Asymmetric Abelian Avalanches and Sandpiles Page 12